if !d.name eq 'PS' then begin
    device,xsize=17,ysize=30,yoffset=3
    !p.thick=1 & !x.thick=1 & !y.thick=1
end

default,png,0
!p.multi=[0,2,2] & !p.charsize=1.3 
!x.margin=[7,0] & !y.margin=[3,3]
; field to be ploted
@eu_setup
field=bx
names=['a) ', 'b) ', 'c) ', 'd) ']
;;;;;;;;;;;;;;;;;;;;;
default,ibeg,0
nt=(size(field))[3]
print,nt
nl=64
year = 365.*24.*60.*60.
time = (nxaver*dt*indgen(nt))/year
good=where(time gt 0)
rg=where(r gt 0.74) 
bmax=(max(abs(field[*,rg,good])))/0.33
bmin=-bmax
lev=grange(-1.0,1.0,nl)
year = 365.*24.*60.*60.
fo = "(F5.1)"
ii=0

; polodial field lines
comp_r = dblarr(l,m)
comp_theta = dblarr(l,m)
lines_r=dblarr(l,m)
lines_theta=dblarr(l,m)
rsinth=fltarr(l,m)
isinth=fltarr(l,m)
r2d=fltarr(l,m)
for k=0,l-1 do r2d[k,*]=rr*r[k]
for k=0, l-1 do begin
  for j=1,m-2 do begin
    isinth[k,j]=1./sin(theta[j])
    rsinth[k,j]=rr*r[k]*sin(theta[j])
  endfor
endfor

;cpos=[0.23,0.08,0.99,0.9]
bpos=[0.20,0.67,0.21,0.83]

 xlength=0.19
 xinterval=0.05
 bini=0.10
 bint=0.24
 xpos1=xinterval
 xpos2=xinterval+xlength
 cpos=[xpos1,0.04,xpos2,0.94]
; bpos=[bini,0.31,bini+0.01,0.63]

for i=235,340,32 do begin
 time = 20.+nxaver*243.*dt/year + nxaver*dt*double(ii*32)/year


 ;; PLOTING
 contour,smooth(transpose(reform(field[*,*,i])),2),x[*,*],y[*,*],$
         /fi,nl=64,/iso,lev=lev,xstyle=4,ystyle=4,$
         title=names[ii]+'!8t=!6'+string(time,FORMAT=fo)+' !6yr'
;         title='!8B!D!7u!N!8 / B!D!6eq!6!N, !8t=!6'+string(time,FORMAT=fo)+' !6yr'
 if (i eq 243) then begin
 colorbar,range=[min(lev),max(lev)],/vertical,xaxis=2,$
         ytickformat='(F6.2)',yticks=4,pos=bpos,charsize=1.5,$
         ytickv=[min(lev),0.5*(min(lev)+max(lev)),max(lev)],$
         title='!8B!D!7u!N!6 (T)'
 endif
; poloidal field lines 
 br=transpose(reform(bz[*,*,i]))
 bth=transpose(reform(by[*,*,i]))
 comp_r= br * rsinth
 comp_theta = bth * r2d
 for j=0,m-1 do begin
  for k=0,l-3 do begin
    kk=indgen(k+2)
    lines_theta[k,j] = int_tabulated(rr*r[kk],comp_theta[kk,j],/double)
  endfor
 endfor
 lines_theta[0,*]=lines_theta[1,*]
 lines_theta[l-2,*]=(4.*lines_theta[l-3,*]-lines_theta[l-4])/3.
 lines_theta[l-1,*]=(4.*lines_theta[l-2,*]-lines_theta[l-3])/3.
  
 for k=0,l-1 do begin
  for j=m-1,2,-1 do begin
   jj=indgen(j)
   lines_r[k,j-1]=int_tabulated(theta[jj],comp_r[k,jj],/double)   
  endfor
 endfor
; lines_r[*,0]=-lines_r[*,1]
; lines_r[*,m-1]=-lines_r[*,m-2]

;  lines = rsinth * 0.5 *( isinth*lines_r - (1/r2d)*lines_theta )
;  lines = lines_r
  lines = -rsinth*(1/r2d)*lines_theta
; end of poloidal field lines
 nls=17
 lev_lines_p = grange(-0.5,max(lines[*,*]),nls)
 lev_lines_n = grange(min(lines[*,*]),0.5,nls)
 contour,smooth(lines,2),x,y,/over,lev=lev_lines_n,$
        max_value=0.,min_value=min(lev_lines_n),c_thick=1,c_linestyle=2,color=cgColor("white")
 contour,smooth(lines,2),x,y,/over,lev=lev_lines_p,$
        max_value=max(lev_lines_p),min_value=0,c_thick=1,c_linestyle=0,color=cgColor("white")
;
 plots,circle(0,0,0.96),thick=3
 plots,circle(0,0,0.61),thick=3
 plots,circle(0,0,0.718),thick=2,li=2
 oplot,[0,0],[-0.96,-0.61],thick=2,li=0
 oplot,[0,0],[0.61,0.96],thick=3,li=0
;; END OF PLOTTING

;xpos1=xpos2+xinterval
;xpos2=xpos1+xlength
;bini=bini+bint
;cpos=[xpos1,0.04,xpos2,0.94]

 print,i
    if (png EQ 1) then begin 
     istr2 = strtrim(string(ii,'(I20.4)'),2) ;(only up to 9999 frames)
     print,"writing png",'./img_'+istr2+'.png'
     image = tvrd(true=1)
     imgname = './img_'+istr2+'.png'
     write_png, imgname, image, red, green, blue
    endif
 wait,0
 ii ++
endfor

!p.multi=0
end

